opts_knit$set(progress = TRUE, verbose = TRUE, root.dir = "~/Documents/git/ash/paper/Rcode")
require(ashr)
## Loading required package: ashr
## Loading required package: Rcpp
require(qvalue)
## Loading required package: qvalue
require(fdrtool)
## Loading required package: fdrtool
require(mixfdr)
## Loading required package: mixfdr
require(locfdr)
## Loading required package: locfdr
## Loading required package: splines
require(ggplot2)
## Loading required package: ggplot2
##
## Attaching package: 'ggplot2'
##
## The following object is masked from 'package:qvalue':
##
## qplot
load("sim1_out.RData")
#look into methods for estimating pi0 conservatively
# par(mfcol=c(3,3))
plot_ecdf = function(sims) {
for (i in 1:length(sims$beta)) {
plot(ecdf(sims$beta[[i]]), xlim = c(-6, 6), main = paste0("iteration ",
i))
x = seq(-6, 6, length = 1000)
lines(cdf.ash(sims$betahat.ash.n[[i]], x), col = 2, lwd = 2)
lines(cdf.ash(sims$betahat.ash.u[[i]], x), col = 3, lwd = 2)
lines(cdf.ash(sims$betahat.ash.true[[i]], x), col = 4, lwd = 2)
}
}
plot_ecdf(simres1)
plot_ecdf(simres2)
Figure to show that estimated betahats are not so different
plot(betahat.ash.u[[1]]$PosteriorMean, betahat.ash.n[[1]]$PosteriorMean)
## Error: object 'betahat.ash.u' not found
abline(a = 0, b = 1, lwd = 2, col = 2)
## Error: plot.new has not been called yet
QUestion: is the null-biased prior maybe a little too conservative? Answer: log likelihoods don't suggest they are
# hh.ashtrue = hh.ashz hh.ashtrue$fitted.g$pi =
# c(2/3,1/15,1/15,1/15,1/15,1/15) hh.ashtrue$fitted.g$mean = c(0,0,0,0,0,0)
# hh.ashtrue$fitted.g$sd = sqrt(c(0,1,0.2,0.4,0.8,3))
# loglik(hh.ashtrue,betahat,sebetahat) loglik(hh.ashz,betahat,sebetahat)